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ABSTRACT 

A numerical study of laminar film boiling by finite- 
difference method has been carried out. The finite-difference 
form of the basic conservation equations and interface condi- 
tions are solved by the false-transient method. The effect 
of radiation and the variation of thermophysical properties 
with temperature are considered in detail. Ratification of 
the present work is done by comparison with earlier studies. 
Velocity j temperature, vapour layer thickness and average 
Nusselt numbers are obtained for water at a pressure of one 
bar for various degrees of superheating and subcooling. The 
effect of radiation on vapour layer thickness and total heat 
transfer rate are illustrated. These results are also compared 
with experimental data on saturated and subcooled film boiling 
available in the literature. 



CHAPTER I 


INTRODUCTION 


1.1 Film Boiling ; 

The name, film boiling, has been given to that type of 
boiling which occurs when a complete vapour film exists between 
the heated surface and the boiling liquid. Heat exchange 
processes accompanied by film boiling are wide spread in such 
fields of technology as,' metallurgy, during thermal treatment 
of metals? radio electronics, in temperature control of 
electronic equipment, and power production etc. Also, recent 
safety measures on nuclear reactors stipulate that the condition 
of film boiling should thoroughly be investigated. 

During the quenching of materials, a very high tempe- 
rature material is suspended inside a liquid. During the 
initial period a complex vapour film exists and the cooling of 
material is by film boiling heat transfer. The rate at which 
the material is cooled initially is important in determining 
its structure and hence the film boiling plays a vital role 
there . 

In radio-electronics film boiling becomes important 
in temperature control of electronic equipment. In cryogenic 
engineering we come across liquids of very low boiling point. 
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Film boiling becomes important if the cryogenic liquids 
come across surfaces of moderate temperatures. Operation of 
jets or rockets frequently involves the contact of a boiling 
liquid with hot surfaces. 

1.2 Variation of Thermophysical Properties with Temperature 
and Pressure ? 

In film boiling the temperature difference across the 
vapour film is inevitably so large that the temperature dependence 
of properties in the vapour film is severe. This situation 
may be easily understood if one thinks about large temperature 
difference between the saturated liquid and the surface and 
the general steep temperature variation of properties of water 
vapour in the proximity of saturation temperature. Also, the 
properties vary with pressure and if the S3?'stem is under pres- 
sure, the heat transfer results might be severely affected 
by the variation of properties. 

The properties that are normally used for calculation 
of heat transfer are calculated either at saturation tempera- 
ture or at some mean temperature between the wall and the 
saturation temperature. Evidently, properties calculated at 
saturation temperature will give much error because of large 
temperature difference between the wall and saturation tem- 
perature. This can be easily seen because the properties at 
the wall temperature may be more than twice as large or small 
than the properties at the saturation temperature. The 
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properties calculated at the film temperature may be fair 
enough if all the properties vary linearly with temperature 
and the temperature profile is linear inside the vapour-layer. 
But the density of vapour-varies inversely with absolute 
temperature (at low pressures) and other properties vary 
according to some power— law and hence the calculation using 
the properties at the film temperature will also be in error. 
The properties vary nonlinear ly with pressure and hence if the 
system is under pressure the constant property model will be 
in error. 

1 . 3 Effec t of Radiation ; 

Under the film boiling conditions, the rate of heat 
transfer from a heated surface to a boiling liquid is effected 
by both convective - conductive transport and by radiative 
transport. The radiative heat transfer becomes increasingly 
more important as the surface temperature increases relative 
to the bulk liquid temperature. 

There are two processes by which the radiative transport 
may affect the film boiling heat transfer. The first is a 
direct transfer between the heated surface and the liquid. 

The second is the absorption and emission of radiation in. the 
vapour film which lies hetween the surface and the liquid. In 
some cases the vapour may have no absorption or emission bands 
which lie in the wavele ' tth range of importance to radiative 
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heat transport. Correspondingly the only radiation effect 
is the direct transfer from the surface to liquid. Typical 
non-participating gases are Oxygen, Hydrogen, Nitrogen and so 
forth. Many vapours including steam have absorption and 
emission bands in the thermal range. But the vapour-film 
thickness is usually very small and we can consider it 
optically thin for radiative transport in most cases. This 
way, we can neglect the effect of the participating gas in 
radiative heat transfer. 

1.4 Review of Literature ; 

Bromley (1950) was the first person to propose a model 
for the problem of film boiling in the presence of radiatively 
non-participating vapour. He employed a simple model which 
neglected inertia forces and superheating within the vapour- 
film and additionally assumed a zero longitudinal vapour 
velocity at the liquid- vapour interface. Also, the properties 
were assumed constant at film temperature. As discussed in 
his paper, Bromley attempted to incorporate the effects of 
surface to liquid radiation into the energy equation for the 
conductive-convective transport. The resulting differential 
equation appeared not to be readily solvable. As an alter- 
native he used qualitative arguments to derive a heat transfer 
coefficient for the combined convection and radiation process. 



P.W. McFadden.and R.J.'Grosh (1961) analysed the prob- 

i 

lem of film belling as ‘a boundary layer problem considering 
the vapour as a compressible fluid. Only the variations of 
density and specific heat were considered with temperature 
and pressure. Other properties were taken constant at film 
temperature. Also, zero-interface velocity was assumed and 
radiation was not dealt with. 

Later J.C.Y. Koh (1962) made an analysis of the problem 
in which the shear stress and vapour-velocity at the vapour- 
liquid interface were taken into account, neglecting the 
variation of properties and radiation part. Sparrow and 
Cess (1962) formulated the problem within the frame work of 
boundary-layer theory in the presence of subcooled liquids, 
with longitudinal velocity at the vapour- liquid interface 
assumed as zero. 

A thorough investigation with regard to radiation was 
made by Sparrow in 1964. Consideration was given both to the 
direct radiation between the heated surface and the liquid and 
to the emission and absorption of radiation in the vapour 
layer. He found that the direct surface-to-liquid radiation 
can appreciably increase the film boiling heat transfer . A 
quantitative criterion was deduced which states the conditions 
under which the effects of surface-to-liquid radiation are 
significant. The results of his analysis (with radiatively 
participating gas) indicated that the effects of a radiatively 
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participating vapour on heat-transfer were negligible. In his 
analysis he assumed zero longitudinal- velocity at the vapour 
liquid inter face. 

K. Nishikawa and T.Ito (1966) solved this problem in 
the presence of subcooled liquids in which shear stress and 
longitudinal velocity at the interface were taken into account. 

In their analysis, properties were assumed constant and radia- 
tion was neglected. 

Recently, Nishikawa et.al. (1976) have investigated 
the variably thermophysical property problem and have showed 
that the constant property problem gives results which differ 
widely from that of variable property problem for large tem- 
perature differences between the wall to liquid and for higher 
pressures. Radiative energy transfer was neglected in the 

analysis. 

Sabhapathy (1980) studied the variable property problem 
with radiation using finite difference method. This is the 
first known study of film boiling, in which the governing 
partial differential equations were solved using finite diffe- 
rence techniques. His analysis can be improved significantly 
because of the following reasons. In that analysis, some signifi 
cant terms which arose due to the transformation of coordinates 
(See Sec. 2.4 and Appendix A in this work) were neglected. 

More importantly , the convergence criterion he adopted was 
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not reliable, when the time steps were very small (See Sec. 3.1 
in this work). Also, he employed the false— transient technique 
to solve the governing time-dependent partial differential 
equations, except the equation governing the boundary layer 
thickness. He solved the last equation using- marching tech- 
niques. With this approach, however, he was forced to fix 
the boundary layer thickness at the first streamwise location 
to control the instabilities. 

1 . 5 Intr oduction to the Present Work ; 

In the present work, the effect of property variation 
and radiation on film boiling are considered together . The 
finite difference form of the basic conservation equations 
and interface conditions are solved by the false— transient 
method. The effect of radiation on heat transfer coefficient 
is clearly brought-out. Heat transfer coefficient is calculated 
for various degrees of superheating and subcooling. The 
accuracy of the present method is checked by comparison 
with results of Nishikawa et.al. (1976) for the case of no 
radiation. The results are also compared with experimental 
data on saturated and subcooled laminar film boiling available 
in the literature. 


f 



CHAPTER II 


ANALYSIS 


In this chapter, the finite difference form of the 
basic conservation equations in laminar film boiling from a 
vertical plate are developed. The symbols used for various 
parameters and subscripts are indicated under nomenclature. 

2.1 Physical Model ; 

The physical model and the coordinate system for the 

plane vertical heated plate are shown In Fig. (2.1). The 

heated plate whose surface is kept at a uniform temperature 

T is submerged vertically in a stagnant boiling liquid whose 
w 

temperature is lower than the saturation temperature T vg by 
the degree of subcooling aT^. 

Important assumptions made for the derivation of the 
fundamental equations of conservation laws are; (l) the 
vapour — liquid interface is smooth and held at saturation 
temperature, (2) Boundary ~ layer approximations are valid 
in the vapour film surrounding the heated plate and the liquid 
adjacent to it, (3) vapour is assumed to be a non-participat- 
ing gas while considering the radiative energy transfer from 
the wall to the interface. 
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By the physical model and the coordinates described 
above the fundamental equations for laminar film boiling for 


vapour and liquid boundary layers are given as follows 
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The first three equations correspond to the vapour 
layer. Equation (2.1) representing the continuity, equation (2.2) 
representing the momentum conservation, and equation (2.3) 
representing the energy conservation in the vapour layer. 

The last three equations represent the continuit}?- , 
momentum conservation and energy conservation of the liquid 


layer respectively. 

These equations must be compatible with the following 
relations at the vapour-liquid interface in steady state. 
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Equation (2.7) represents the conservation of mass at the 
vapour-liquid interface. The continuity of tangential velocity 
at the interface is represented by equation (2.8). It is 
nothing but the no slip condition at the interface. Equation 
(2.9) represents the matching of shear stress at the interface. 
Equation (2.10) represents the conservation of energy at the 
interface. The continuity of temperature is represented by 
equation (2.11). 


Other conventional boundary conditions can be written as, 

at y = 0 U v = V v = 0 and T y = T w (2.12) 

as y - co U L - 0 and T L - T^ (2.13) 

Equation (2.12) represents the no slip condition at the wall 
and the temperature specification. At large distance from the 
wall the liquid is undisturbed by the presence of the wall and 
this is represented by equation (2.13). Energy is transferred 
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i 

from the wall to the liquid by conduction - convection and 

i 

radiation. The term q R represents the radiation heat-transfer 
term. qp has been evaluated assuming the wall is not black 
and its emissivity is s w and interface is assumed to be black. 

2.3 Selection of Numerical Technique : 

This type of problem is usually solved by similarity 
transformation. In the presence of radiation interaction 
similarity transformation is not possible. So s the partial 
differential equations are solved numerically by using finite- 
difference techniques. 

These equations can be solved either by marching 
technique (since the equations are parabolic in x) or by a 
time— dependent approach such as the false— transient technique. 
An early work in which both the stead}^ state and the time- 
dependent approaches were used was that of Hung and Macagno 
(Roache, 1976). They found that the time-dependent equations 
were easier to handle and more stable , a conclusion confirmed 
bv many workers since then. Eventhough the philosophy of the 
time-dependent approach is attractive * it is natural to ask 
why should one bother with the transient solution in those 
cases where the only interest is in the eventual steady state 
solution. Why not set 3w/ 3t = 0 and just work with the steady 

flow equations ? The conclusion obviously depends on the 
simplicity of the time-dependent method used. For reasonably 
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simple time-dependent methods, the flexibility of being able 
to achieve the transient solution if desired, is also attrac- 
tive. More important, the time— dependent approach does not 
presume the existence of a steady state solution, which indeed 
may not exist. It also gives more physical insight into the 
nature of the equations. This is important in a two-phase 
problem which has complexities on account of interaction 
between liquid and vapour. 

An attempt was made by Sabhapathy (1980) to solve these 
equations by marching technique but he was unaQle to control 
the instabilities even with very small mesh sizes. We solved 
the equations using false-transient technique. Here the 
unsteady state problem is marched in time till a steady state 
solution is achieved. 

2 . 4 Selection of Grid System : 

There are two ways of introducing grid points into the 
coordinate system. In the first method as shown in Fig. (2.2*1), 
the x— y coordinate system can be used. If the numerical 
solution is sought in this x-y coordinate system then it has 
one important defect? the coordinate system (Cartesian) will 
not, in general, conform to the shape of the boundary-layer. 

As the boundary layer grows, more grid points must be added or 
be carried from the outset of the problem. Thus the scheme 
involves testing for the outer edge of boundary layer, at 
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every downstream step hy checking the change in the value of 

V/. . at each cross— stream step? the outer edge being defined 
il 

by the point where j (-~y ) a I is less than a prescribed 
tolerance. Such a method is uneconomic of computer storage 
and time. Also, there are grids of curved shape near the 
vapour-liquid interface, which need special equations to be 
written at these points. A better technique is to emploir a 
coordinate system that grows with the boundary layer. 

In the second method, as shown in Fig. (2.2»2), instead 
of using y as the cross— stream coordinate, we could use the 
normalized distance n = y/fi > <5 being the vapour boundary 
layer thickness at any x, such that n = y/$ = 1 corresponds 

to the edge of the vapour boundary layer. In this way, the 
finite difference grid for 0 <_ y/ 6 <_ 1 would always fit the 
boundary layer region and the grid points selected initially, 
at x = 0, say, would lie within the region of interest for 
all values of x. As the vapour liquid interface always made 
to correspond to a grid point there can be no problem of 
curved grids. In the present work, numerical solution is 
obtained in x - n coordinate system. The following section 
deals with the changes in the equations of conservation due to 
transformation of coordinates from x - y to x - n system and 
the addition of transient term to proceed with the false 


transient scheme. 
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2.5 


the x - n Coordinate System : 


'State Problem in 


The basic conservation equations- given in the x - y 
coordinate system in Sec. (2.2) are to be modified in the 
x — -n coordinate system, because n - y/s is a function of 
. both x and y. The equivalent forms of first and second par- 
tial derivatives of a variable in the new (x - n) coordinate 
system are given in Appendix A. The governing equations of 
the unsteady state problem that are solved are as follows. 
For the vapour side, 
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The compatibility conditions at the interface ares 
At n = Is 
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The mass conservation equation at the interface can be derived 
(Appendix B) and is given as follows. 
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Similarly the interface energy equation in unsteady state can 
be stated as follows (See Appendix B for derivation). 
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The conventional boundary conditions are; 
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as 


n = 0 O v - V- 0 


and T v = T w 
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2 . 6 Dimensionless form of the -Equations ;; 

The governing equations were non-dimensionalizea and 
the dimensionless variables are given as follows. 
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The non-dimensional numbers evolved from this non-dimensiona- 
lisation are given as follows. 
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Peclet number of vapour 

Jacob number of liquid 
Peclet number of liquid 
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From now on in the dimensionless form of the governing 
equations primes will be removed for brevity. From these 
final governing equations the finite difference equations will 

be developed. 

Momentum Equation on Vapour Side ; 
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Liquid Side 


Momentum 
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Che compatibility^ conditions at the interface arei 


at n = 1 
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Mass f!nnsprvation at the Interfaces 
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Interface Energy Balance Equation ; . 
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2.7 Numerical Details and Method of Solution s 

The length of the plate over which the numerical solu- 
tion obtained is 0.015 m, because flow is generally not laminar 
above this height. Although this may seem like an unrealisti- 
cally small height, the results obtained here can be compared 
with experimental results on film boiling from thin wires 
( < 10 mm dia.) where the curvature effects can be neglected. 
Thus, the number of grid points K chosen along the x-direction 
deter mi nes the space step in the x-direction ( * 0.015/K). In 
the present work K was chosen to be 15 « 

The number of grid points N, chosen along cross-stream 
direction in the vapour boundary layer was 10, the N-th grid 
point corresponding to the vapour— liquid interface. 1=0 
corresponds to x = 0, J = 0 corresponds to the wall and M 
corresponds to the liquid at <». 

The liquid boundary layer is a conventional boundary 
layer, its thickness can be defined by arbitrarily defining 
the edge of the same. The number of grid points was chosen 
to be 30 ( * M - N). This was found to be sufficient - 
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and was also verified from the analytical work of Frederking 

* 

et.al. (1964) - 'See Fig. 2.5. 

Method of Solution ? 

The momentum equation and the energy equation, given 
above for both gas and liquid, are parabolic in time, contain- 
two independent space coordinates and are coupled through the 
buoyancy term. These equations along with the continuity 
equation and the interface compatibility conditions are 
solved, with, the corresponding boundary conditions and initial 
conditions, by finite difference techniques, employing an 
explicit, two level method. It is explicit because all the 
quantities needed to calculate values at the new time step by 
the finite difference equations are known. It is two-time level 
because only two-time levels are involved in the calculation. 

The partial derivatives in partial differential equa- 
tions are approximated by suitable finite difference expres- 
sions (central difference scheme for diffusion terms and 
upwind difference scheme for convection terms). This proce- 
dure leads to a set of algebraic equations from which the 
solution for the present time step can be obtained. 

By making the grid spacing s sufficiently small, the 
solution obtained by this procedure, is a sufficiently close 
approximation to the exact solution, if the scheme is conver- 
gent, stable and consistent. 
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1 The sequence for advancement from an old time level 
t to a new time level t + At is given as follows. 

(1) Velocity (U ) is calculated in the vapour region for 
J = 1 to N - 1 from the finite difference form of the vapour 
momentum equation (2.27). 

(2) Velocity (U^) is calculated in the liquid region for 
J = N+l to M-l from the finite difference form of the liquid 
momentum equation (2.30). 

(3) Velocity (U = U, ) is calculated from the finite 

v int ^int 

difference form of interface shear stress condition (2.34) « 

(4) Temperature £T ) is calculated in the vapour region 
for J = 1 to N - 1 from the finite-difference form of energjr 
equation (2.28) . 

(5) Temperature (T^) in the liquid is calculated in the 
liquid region for J = N+l to M-l from the finite difference 
form of the energy equation (2.31). 

(6) Velocity (V ) is calculated in the. vapour region for 
J = 1 to N from vapour continuity equation (2.29). 

(7) Boundary-layer thickness ( 3 ) is calculated from the 
finite difference form of the interface energy conservation 
equation (2.36). 

(8) Velocity (V T . ) is calculated from the finite difference 

int 

form of the interface mass conservation equation (2.42). 

(9) Velocity (V^) is calculated in the liquid region from 
the finite difference form of liquid continuity equation (2.39) - 
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Finite Difference Equations ; 

FDE’s were developed by difference approximation to the 
partial derivatives of the governing equations. These FDE's were 
used to find variables like velocities s temperatures and the 
vapour boundary layer thickness. 


Finite Difference Approximation to Calculate the Velocity (U v ) 
from the Equation (2.27) for J = 1, N-ls 
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Finite Difference Approximation to Calculate the Velocity (U L ) 
from the equation (2.30) for J = N+l to M— 1. (Notes In the terms 
appearing with * ’s index 3 is to be replaced by j -1 if the 
velocity V L is positive at that location-upwind scheme). 
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Finite difference approximation to calculate the interface 
tangential velocity (U = U, ) from the equation (2.34-). 
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Finite difference approximation to calculate the temperature 
(T^) from the equation (2.38) for J = 1 to N— Is 
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Finite difference approximation to calculate the temperature 
(T^) from the equation (2.31) for J = N+l to M-l (Notet In terms 
appearing with replace j by j-1 if is positive). 
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Finite difference approximation to calculate the velocity (V v )' 
from the equation (2.29) for j = 1 to N: 



Finite difference approximation to calculate the boundary 
layer thickness ( ) from the equation (2.36). 


(2,42) 
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Finite difference approximation to calculate the velocity (V^) 
at the interface from the equation (2.35).' 
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Finite difference approximation to calculate the velocity (V L ) 
in the liquid region £rom the equation (2.32) for J * N+l to 
M— It 
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It is worth noting that in the above FDE's, all the terms 


V , K » C , P T , v T , K t and C representing thermophysical 
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properties at any location (I,' J) are also dimensionless 
variables. Functional variation of these properties with 
temperature is given in > Appendix C. 

The following information might be useful when iterat- 
ing the finite difference scheme given above. It was mentioned 
earlier that the number of grid points' in ; x-direction was 15. * 
When the FDE's were iterated at all' the 15 grid points, with 
the initial conditions, the following observation was made. 

All the variables at x - locations nearer to the leading edge 
reached steady state .in less number of iterations when com- 
pared to the variables at far downstream. Variables nearest 
to the leading edge reached steady state in smallest number 
of iterations. The reason for this can be drawn as follows. 

For a variable at any x-location, its domain of 
dependence includes a part of upstream region and not down- 
stream region, which is a characteristic of parabolic flows. 
Since the variables at the leading edge are fixed, the vari- 
ables nearest to the leading edge reach steady state in 
smallest number of iterations, because significant number of 
variables in its domain of dependence were fixed. Based on 
this idea, variables at the first x-location (correspond to 
1=1) were iterated till steady state was reached making use 
of variables at leading edge (I = 0). After all variables at 
1=1 reached steady state, variables at next x-location (I = 2) 
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were iterated making use of variables at previous x-location 
(I = l)- Thus this procedure was repeated upto the last 
x-location. 

The main advantages in this method of solving the 
problem are c 

i) Reduction in CPU time, because at each x-location 
variables converged to steady state in same number of itera- 
tions (~ 2000) . 

ii) Computer storage is reduced by a huge amount. In this 
method of solution only the following variables are to be 
stored? 

a) Variables at the current x-location (say I) at 
which the variables are being iterated. 

b) Variables at the previous x-location (1-1). This 
way of obtaining the solution reduced the run t im e 
charges substantially. The following schematic 
diagram might help to quickly understand the above 
mentioned method of solution. 

Choice of Initial Profiles ? 

Vapour velocity in the x-direction (U ) and the vapour 
boundary layer thickness (<$) were substituted from Bromley's 
analysis (1950) at all values of x. 




Schematic Diagram 
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(See Sec. 3.2 for further details about the effect of each of 
these profiles on stability and convergence of the method) 
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Vapour velocity in the y-direction was integrated from the 


continuity equation and the U velocity profile (2.47). 
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CHAPTER III 


NUMERICAL EXPERIMENTS 


This chapter is concerned with the difficulties encoun- 
tered in obtaining numerical solution to the set of differen- 
tial equations presented in the last chapter. The information 
about numerical stability and convergence of the scheme can 
not be easily understood from theoretical analysis because 
of the complicated set of equations with their coupled nature. 
The stability and convergence aspects can, however, be studied 
from numerical experiments, the discussion of which is the 
subject matter of this chapter. Basic problems of convergence 
to a steady state solution and rate of convergence are dis- 
cussed at length. 

3 . 1 Convergence Criteria ; 

The word convergence is used in two different senses, 
iteration convergence and truncation convergence. The latter 
is well known to us, which has to do with the convergence of 
the solution of FDE’s to the solution of PDE's as all step 
sizes tend to zero. Iteration convergence refers to finally 
arriving at a steady state solution, of a set of PDE’s via 
iteration. 
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There are several problem dependent firry, suggestions 
which we can make about attempting to establish iteration 
convergence. Out of those, two widely used criteria are 
discussed here. 

a) Relative error criterion is usually stated as follows: 


max 

il 



< e | ^ ° 


(3.1) 


Here W 1 ?. can be any variable like velocity, temperature or 
boundary layer thickness at the space point 1,3 and at time n. 
This criterion has some weakness in explicit transient methods 9 


when time step is very small* In many transient schemes $ time 
step is governed by stability requirements* If the time step 
turns out to b every small, then one has to think of the 
numerator value in the above criterion (3«l)* I n such, situa- 


tions, we can take the time step almost equal to zero* Conse 
quently ¥ might not be significantly different from 
and yet making the whole term on the LHS less than the pres- 
cribed value e, even in very, initial stages of iterative scheme. 
This might certainly lead us to a premature termination. To 
avoid this, we should constrain e to take very small values, 
not greater than 10 ^ or go for some other criterion. 


Thus the above criterion might fail to lead us to 
steady state in transient schemes when time step is small .When 


4, I * . 

QW%l «n ; 
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the question of what is real steady state ? is to "be answered, 
the following criterion might answer in a better way. 

b) The other .criterion can be stated as follows 

- W n . MOST SIGNIFICANT 

max | - r i <_ e { term in that } (3»2) 

ij equation 


In the present work e is taken as 0.01 which was found to be 
accurate enough. The word MOST SIGNIFICANT is important 
because, in general all the terms in an equation need not 
be of the same order of magnitude (ex: convection terms in 
momentum equation of vapour are very small when compared to 
viscous and buoyancy terms) . Thus comparing transient term 
with a less significant term in that equation might be a strin- 
gent requirement, if not difficult to achieve in practice. 

3 . 2 Choice of Initial Conditions ; 

To solve the FDB's derived in the last chapter the 
following initial conditions were used for the saturated film 
boiling with no radiation. 

1) A quadratic velocity profile with no shear stress at 
the interface was chosen as U velocity profile in the 
vapour layer. 

i.e. U ~ (2 - ^ (3.3) 

v 6 6 

2) Linear temperature distribution in vapour layer. 
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3) 3) V velocity was obtained by integrating the conti- 
nuity equation on vapour side. 

4) U velocity on the liquid side was a quadratic profile. 


5) 


= U , (1 

lq vapour at 

interface 

6 values were obtained from Bromley's analysis (1950) • 
When FDS's were solved with these initial conditions, 
a steady state solution was not reached. The reason for 
instability was traced to the choice of zero shear 
stress velocity profile in vapour layer. We found that 
the evolution of U velocity profiles from the initial 
ones was strongly depending upon the choice of initial 
velocity profile. Fig. (3.1) shows the possible choices 
of initial profiles and their possible directions of 
evolution, when the equations were solved with zero- 
shear velocity profile, the interface shear condition 
satisfied in the way shown in Fig. (3.1a). Also, the 
solution was not converging to steady state. Then the 
initial U velocity profile was suspected to be an 
inappropriate one and Bromley's zero interface velocity 
profile was substituted as initial profile, 



(3.4) 

L 


U. 


y/s - y/ 6) 


(3.5) 


vapour 

With the above velocity profile, the solution to FDE's 
converged to steady state in 22000 iterations (approxi- 
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mately 20 minutes of CPU time). This is certainlv a slow 

J 

4 

convergence-. 

3.3 The Effect of Initial Value of Boundary-Layer Thickness ; 

It was found that the steady state value of 5 was very 
much different from its initial value obtained from Bromley's 
analysis (See Fig. 3.2). when more refined value of 6 was sub- 
stituted as initial data along with other initial conditions 
for solving subsequent problems (with varying degrees of 
superheating and subcooling) , a further delay in convergence 
was observed which can be explained as follows. 

The delay in convergence was due to slowly developing 
(increasing from zero .everywhere) U velocity profile on the 
liquid side. If large momentum was supplied to the liquid 
from the interface then the velocity profile might quickly 
develop to its steady state value. Hence the V velocity in 
the liquid should be positive. This implies that 6 should 
increase from the initial value because, d^dt represents 
physically the interface velocity (~ V velocity of liquid 
at interface in transient state). But as already mentioned, 
all variables (including . 6 ) , except U velocity in the liquid 
region, reached steady state after 2000 iterations. It 
implies that d£j/dt is almost equal to zero after two thousand 
iterations. And the only way to make the term (dd/dt) positive 
at this stage was to reduce intentionally from its steady 
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value. In this way 6 was perturbed and convergence was 
established in saturated film boiling in 18000 iterations. 

Yet this is not a significant achievement. The slow •onver- 
gence might be due to the use of same non-dimensional time- 
step in both gas and liquid regions. In fact, the gas region 
reached steady state in 2000 iterations and the rest of the 
time was taken up by the liquid. This can be explained as 
follows. 

The time step chosen was because of the Stability 
requirements on gas side. Also, the allowable time step was 
found to be a strong function of the degree of superheating 
and subcooling. For decreasing superheating and increasing 
subcooling, the allowable time step was found to be decreas- 
ing (See Fig. 3«3h). And this could certainly increase th.e 
number of iterations needed for reaching steady state solution 
in case of strongly subcooled liquids. But it was noticed 
that, as we proceed further seeking solution at downstream 
regions, time step can be increased gradually by a signifi- 
cant amount (at least by 4 or 5 times the value of At at 
leading edge). See Appendix D for a brief explanation. This 
can help us in reducing the number of iterations atleast at 
downstream regions. 

The allowable time steps for varying ax values, which 
can give numerically stable scheme are shown in Fig. (3.3a). 
The effect of Ay on time step is very important when compared 
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to the effect of ax. The order of magnitude of each effect 
on time step is given in Appendix D, 

No instabilities were observed on the liquid side 
at this time-step (See Appendix D for a brief discussion). 

Thus the time step used on the liquid earlier, which was 
desired by stability on gas side, was too small for liquid 

m - 

and it could not come-up to steady state in 2000 iterations. 
Based on this, one may conclude that it is desirable to have 
two different time steps in order to take advantage of the 
higher At that can be used in the liquid region. This method 
will create difficulties with regard to satisfaction of 
interface conditions. The other possibility is to decouple 
gas and liquid equations and run them separately and couple 
them later. This approach is discussed next. 

As the two regions are coupled by the interface condi- 
tions, one has to fix a set of values (for velocities) at 
the interface, before decoupling can be done. These values 
were taken as the values obtained in the last iteration when 
both gas and liquid regions were coupled. Thus the two regions 
were, decoupled and a relatively large time step, approxi- 
mately 100 times the allowable time step on the gas side, was 
used on the liquid side. Convergence was established on the 
liquid side only in 100 iterations, with this large time step. 
The advantage of this decoupling is that we can force the slowly 
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converging liquid region to arrive at steady state solution 
separately in a small number of iterations. But this steady 
state on the liquid side may not satisfy the interface 
conditions. To achieve this, the two regions were coupled 
for 1000 more iterations and steady state was ex amin ed at 
each and every point. At this stage, all the variables did 
not satisfy the convergence criterion. Again the liquid region 
was decoupled and the same procedure was repeated till conver- 
gence is established at all points when both gas and liquid 
regions were coupled .With this approach we arrived at steady 
state solution in four decouplings (see Schematic diagram 3.4). 

A more significant achievement was possible with the 
use of an over-relaxation technique for the liquid region. 

The following equation shows the use of relaxation parameter 
in the liquid momentum equation. 

U n+1 - u n +0 At T coir tribution due to convection, -i 

* L buoyancy and diffusion terms J 

( 3 . 6 ) 

0 is the relaxation parameter. When the value 0 was chosen 
to be in the range of 100 - 400 (over-relaxation) it was 
found that the FDE's converged to steady state within 2000 
iterations (without any need for decoupling liquid and vapour). 
This is a major improvement as compared to 18000 - 22,000 ite- 
rations required by previous techniques. If a value above 400 
is used, then the equations become unstable .With this 



Are the interface condi- 
tions satisfied 


Schematic Diagram 3.4 
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accelerated convergence the CPU time is 4 minutes for 15 grid 
points in the x-direction. 

3.4 Subcooled Film Boiling ; 

The numerical solution of conservation equations for 
the case of subcooled film boiling is somewhat more difficult 
than that for saturated film boiling because of the presence 
of energy equation on the liquid side in the former. In the 
case of subcooled film boiling, the liquid at the interface 
removes some energy by conduction and convection for heating 
up the liquid thermal boundary layer. Hence less energy is 
available for evaporation at the interface. This leads to 
a reduction in vapour boundary ' layer thickness. In a numeri- 
cal scheme, if this reduction occurs rapidly then it can lead 
to instability as explained below.. 

Velocity V^in the liquid at the interface is calculated 
from mas : s conservation equation (2.35). Hence, this velocity 
is dependent on (l) rate of change Of boundary . layer thickness 
(d6./dt), (-2) shape of the vapour boundary layer (d5/dx), 

(3) rate of liquid disappearance because of evaporation.' 

At low pressures the density of liquid is much larger than that 
of vapour and hence the third term above makes negligible 
contribution. In the numerical scheme the initial velocity Vy 
on the liquid side was based on the results of saturated 
film boiling. If the degree of subcooling is 1°K then the 
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finite difference equations converged satis f act or ily. In 
tne initial stages, the time step was chosen to be half that 
of saturated film boiling. At higher degrees of subcooling 
(around 5°K) the numerical scheme showed tendenc 3 r to become 
unstable. This can be explained as follows. 

In the initial stages the velocity V is deter- 

j 

mined primarily by the rate of change of vapour boundary- 
layer thickness with time (d6/dt). In the case of 1°K sub- 
cooling the velocity V y was small (~ 0.3 m/S) because 
( ds/ dt) was small. On the other hand in the case of 5°K 
subcooling the velocity V v became very large (~ I400m/S) 
and this resulted in instability. When the velocities in 
the liquid region are ver 3 ^ large, the present anal 3 ^sis will 
not be valid (for instance, bounda^r layer approximations 
will not hold good when the cross-stream velocity is 
comparable to streamwise velocity). 

One of the reasons for the large d<$/dt in higher 
degrees of subcooling was traced to a significant change in 
temperature in the liquid. In order to bring only a small 
change in temperature with time an under relaxation technique 
was employed to energy equation which was similar to the one 
described in Secs. 3. 3 and 3.5, but with a value of 0 = 0.01 
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(under relaxation). But this became not successful. Hence? 
it was decided to neglect the term dfi/dt totally in order to 
® liroinate this instability. Nevertheless the set of 
equations are still quite valid in steady state because 
db/dt is zero in steady state. This was even verified from 
the saturated film boiling results. With this modification 
no further instability was observed and the iterated solution 
converged to steady state. 

3 • 5 Over-Relaxation Technique ; 

While solving the FDE’s for the case with subcooling 
an over-relaxation was used which accelerated the convergence 
on the liquid side. The following equation shows the use 
of relaxation parameter in the liquid momentum and energy 
equations. 

wJJ 1 = W?. +0 At [A + B + C ] 

where can be either temperature or velocity on the 

liquid side, and 

^ ~ Contribution due to convection terms In momentum and 

td n j. J _. } energy equations 

B = Contribution due to diffusion terms 

C = Contribution due to Buoyancy term in momentum eon. 

=0 in energy equation 


(See Eqns, 2.30 and 2.31) 
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The value of 0 was varied from 100 to 500 (over relaxation) 
and it was found that for values greater than 400 the scheme 
was unstable. In the present work its value was taken to be 
250. Even though this technique was first developed in 
solving subcooled film boiling, the technique was also employed 
in solvln S FDE's in saturated film boiling for varying degrees 
of super heating and already mentioned in Sec. 3.3. With 
this accelerated convergence the iterated solution reached 
steady state in 2000 iterations with a CPU time of 7.5 minutes 
for fifteen grid points in the x-direction. Also, it was 
found that a significant amount of CPU time was taken up by 
calculation of properties (around 3 minutes of CPU time). 



CHAPTER IV 


RESULTS AND DISCUSSIONS 


The present finite-difference scheme is ratified by- 
comparing heat— transfer coefficient* Velocity and temperature 
profiles (for NO RADIATION CASE) , with those obtained by 
Nishikawa et.al. (1970) who used similarity transformation 
and solved a set of ordinary differential equations. 

In order to examine the effect of variable properties 
and radiation, the degree of superheating, subcooling are 
varied for water at 1 bar pressure. All the relevant proper- 
ties in the governing equations are regarded as temperature- 
dependent when solutions with variable property are to be 
obtained, while in calculation with constant property these 
are evaluated at the saturation or the film temperature (arith- 
metic mean of the wall and saturation temperature). Solutions 
are obtained for various degrees of superheating and subcooling, 
when radiative effects are neglected. By a comparison of the 
results obtained without radiation to those obtained with 
radiation, the effect of radiation on heat-transfer coefficient 
is clearly br ought-out. 
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4.1 Ratification of the Present Work ; 

. The comparison of the values of hx 1 ^ in Fig. (4.1) 
for various degrees of superheating- and subcooling indicates 
very good agreement with those obtained by Nishikawa et.al. (1976) . 
The comparison was done with the values of hx 1 ^ because, 

Nishikawa et.al. showed that hx^^ has the same value for given 

macro-parameters (such as pressure, degree of superheating and 
degree of subcooling) and do not depend on the height x. In 
Fig. (4.2) velocity and, temperature profiles in the vapour 
are compared for a typical degree of superheating ( = 449°K 
and aT-^= 40°K) . The comparison indicates a close agreement 
near the wall and some disparity? - away from the wall in both 
velocity and temperature profiles in the vapour region. This 
disparity away from the wall can partially be attributed to 

the error introduced in the transformation of non-linear 

abscissa - used by Nishikawa et.al. to a more convenient 
t) = y/6 abscissa used in the present work (See Appendix E 
for further details). The above comparisons indicate that the 
present finite difference scheme with grid spacing (ax = 1 mm ■ 

Ay = 0.1 6 ) gives acceptable accuracy. 

^*2 Discussion on V Velocity in Vapour and Liquid Regions ; 

In Fig. (4.3), we see the variation of V velocity with 
distance away from the wall in both vapour and liquid regions. 

As we see in the figure, the velocity is negative all the way 
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in the vapour region and is maximum near the inter face. A 
discontinuity in V velocity exists at the inter face which 
can be explained as follows. 

V velocity in liquid region at the interface is found 
from inter face mass conservation equation (2.35). In steady- 
state, the V velocity is positive near the inter face as 
shown in Fig. (4.3). At the first instant, one might be skepti- 
cal about the validity of this result, because a negative V 
velocity is required to supply liquid at the inter face. But 
the fact that even this positive V velocity can supply liquid 
at the interface is explained based on the curvature of the 
inter face. By considering the streamline drawn in Fig. (4.4a), 
we can conclude the following, (l) Liquid is moving towards 
interface all the way. Thus, it supplies required mass at the 
interface. (2) In the coordinate system shown in Fig. (4.4b) 
velocity of liquid near the interface is positive. 

So, this positive V velocity near the interface is a 
feasible solution. With this, the mass conservation is satis- 
fied as shown in Fig. (4,4c). Other possible way in which 
mass conservation can be satisfied is shown in Fig. (4.4d). If 
this is the way in which mass is conserved, then it says that 
liquid with its density 1600 times the density of vapour (at 
low pressures) is entering vapour-layer in BOTH x and y direc- 
tions ONLY to supply the net amount of vapour leaving that 
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particular differential element, which might not be' physically 
possible unless the order of vapour velocity (or density) is 
very high. 

Thus, the mass conservation with positive V velocity is 
seemingly more correct at low pressures. At this juncture, we 
would like to discuss an important consequence of this new 
result. 

Interface with no curvature (Flat Interface )t For this case, 
d6 is zero at all x. If liquid is to be supplied at the 
interface, then V velocity in the liquid region should be 
negative for this flat interface condition as shown in 
Fig. (4.4e). For the boundary-layer problem d&Aix (one measure 
of curvature of the interface) varies with x, decreases with 
increase in x. With increase in x, the. curvature of the inter- 
face approaches that of the hypothetical flat interface in 
which case, one can expect more region near the interface on 
the liquid side with negative V velocity. This important 
feature was readier observed in the present work. At the 
leading edge (where d6/dx is maximum) in a comparatively large 
liquid region ( 5 points away from interface) V velocity is 
positive and at downstream stations V is positive in a 
relatively smaller region (on^ 3 points away from the interface), 
see Fig. (4.4f). Thus, it shows clearly that we can attribute 
the positive nature of V velocity near the interface to the 
curvature of the interface. 
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4.3 Effect of the Variatibn cf Properties with Temperatures ? 

In Fig. (4.5)? we see the average Nusselt number 

i 

(defined in Appendix F) for different degrees of superheating 
in saturated film boiling. The results of numerical solution 
for variable properties are compared with those for constant 
properties. In the case of constant property solution there 
are two approaches. One is the evaluation of all properties 
at saturation temperature and the other approach is to eva- 
luate them at the film temperature. This figure shows that 
the evaluation of properties at film temperature is accurate 
enough for most applications because, this approach predicts 
heat-transfer coefficients very close to those with variable 
properties at low pressures. 

But it clearly shows that the evaluation of proper- 
ties at saturation temperature could not predict the heat- 
transfer coefficients accurately even at low pressure (l bar). 
Even though, one might always prefer to evaluate properties 
at film temperature (if the mean temperature can be calculated), 
evaluation of properties at the saturation temperature is an 
alternative to start with when the mean temperature can not be 
found out before hand (in case of constant heat flux boundary 
condition at the wall). So, this figure infers an important 
fact that film temperature criterion is preferable to the 
saturation temperature. And the saturation temperature approach 
predicts a lower average Nusselt number at least by 35 percent 
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at higher degrees of superheating when compared with variable 
property approach. 

In Fig. (4.6), we see the. effect of individual property 
variation on heat transfer for saturated film boiling. From the 
figure it is obvious that the effect of variation of specific 
heat is almost negligible at low pressures. All other proper- 
ties that are assumed constant are evaluated at saturation 
temperature. Curve 1 shows the effect of variable viscosity. 
Its effect is to decrease heat transfer and the average Nusselt 
number is in large error at higher degrees of superheating. The 
effect of variable density along with variable viscosity is to 
further decrease the heat transfer as shown in curve 2. There- 
fore, the effect of promotion of heat transfer caused by the 
increase of body force due to the reduced density cannot 
over come the effect of deterioration of heat-transfer caused 
b. 3 ^ the increase of thickness of vapour film. Curve 3, in one 
wa y, shows the effect of variation of specific heat and in other 
way, shows the effect of variation of thermal conductivity, 

Thus, the variation of thermal conductivity, viscosity and 
density could alone predict the heat transfer results accura- 
tely for engineering analysis. In other words, we can neglect 
the variation of specific heat in solving the variable property 
film boiling problem. 
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4.4 Effect of Radiation on Heat-Transfer Results : 

Figure (4.7) shows the effect of radiation on the 
vapour boundary layer thickness for varying degrees of 
superheating and subcooling, In Fig. (4.7a) we see the" effect 
of radiation on the vapour-boundary layer thickness for 
saturated film boiling. Radiation increases the thickness 
of vapour film everywhere and the effect is seen to be more 
predominant as the vapour layer thickness increases. This is 
because, as the vapour layer thickness increases, the convective 
heat transfer goes down and the radiative heat transfer (which 
is constant) plays a relatively more important role. The same 
figure shows that an increase in the degree of superheating 
increases the effect of radiation on boundary layer thickness 
which is easy to understand* Fig. -(4.7b) shows the effect of 
radiation in subcooled boiling. It shows that the effect of 
radiation is small with increasing subcpolfng and decreasing 
superheating. •' 

Fig. (4.8) shows the effect of radiation on average 
Nusselt number (defined in Appendix F) for varying degrees of 
superheating and subcooling. Many important points can be 
inferred from this figure. First of all, it can be seen that 
the effect of radiation is to enhance the total heat transfer. 
This conclusion is consistent with the results predicted by 
Sharrow (1964). Without radiation, Nu decreases with 
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increase in degree of superheating and increases with increase 
in degree of subcooling. It is interesting to note that when 
radiation is included, the Nu increases with degree of 

superheating for saturated film boiling. In subcooled film 
Codling the effect of radiation is less and the average Nusselt 
number decreases initially and increases after some degree of 
superheating. This trend is very much supported by experimental 
data shown in Figs. (4.9a) and (4.9c). 

Radiation affects the heat transfer in the vapour 
layer in two ways. Radiation directly increases the heat- 
transfer by providing an alternative parallel path. It also 
indirectly decreases the heat-transfer by . convection by 
increasing the thickness of the boundary layer. The combined 
effect is to increase the heat transfer as shown in Fig. (4.8). 

Thus, even though, Nu ayerage without radiation decreases 
all the way with increase in degree of superheating, the inclu- 
sion of radiation shows the complicated nature of variation of 
Nu average for var y in g degrees of superheating and subcooling. 

An attempt was made to compare the results obtained 
from the present work for saturated boiling (with radiation) 
with those obtained by Sparrow (1964). Unfortunately, the 
length of the plate considered by Sparrow was 10 cm which 
is unrealistically high for laminar film boiling. However, the 
effect of radiation in terms of the ratio of q/q Q (q is the 
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local heat transfer with radiation and q Q is the local heat 
transfer without radiation) is compared for a typical value of 
x 1*1 cms). Nevertheless, the comparison indicates that 
the trend of the effect of radiation on heat transfer is same 
in both cases, but the above mentioned ratio turned out to be 
1.66 in the present work and sparrow obtained it as 1.92 for 
a degree of superheating equal to 800 K. The assumptions made 
by Sparrow (Sec. 1.4) might be the main reason for this dis- 
parity (See Appendix F for details). 

Bromley in his work (1950), suggested a simple method 
for including the effect of radiation. The validity of that 
criterion is verified with the results obtained in a more 
precise manner in this work (See Appendix F for the verifica- 
tion of Bromle^r's rule for accounting radiation). From the 
values in the table in Appendix F, it can be concluded that 
Bromley s rule for accounting radiation is quite good and 
the disparity is not more than three percent. 

4.5 Comparison with Available Data ; 

Comparison between the present numerical results (with 
full variable property and with radiation) and the experimental 
data (obtained by Nishikawa et al (1976)) is shown in Fig. (4.9). 

The details regarding the transformation of flat plate 
numerical results to horizontal cylinders are shown in Appen- 


dix F. 
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Experimental data plotted in Pig. (4.9a) are for hori- 
zontal cylinder of 13 mm diameter. Theoretical results obtained 
by Nishikawa et.al. (1976) are also plotted in the same figure. 
As it is clear from the figure, data of saturation film 
boiling for cylinder with 13 mm is in good agreement with the 
present numerical results with radiation. At the same time, 
the theoretical results obtained by Nishikawa et.al. (who 
neglected radiation) predicts a smaller average heat-transfer 
coefficient. Appendix F contains the quantitative comparison 
of both these results with experimental data. 

Fig. (4.9b) shows the experimental data for horizon- 
tal cylinder of 6 mm diameter. The results obtained in the 
present work are in small disagreement with the experimental 
data. As the present analysis is done with more sophistication, 
the disparity can be attributed to the curvature effect because 
of the small diameter of the cylinder, since in this case, the 
thickness of the vapour film is comparable to radius of the 
cylinder. 

Experimental data for subcooled film boiling (aT l = 20K) 
are plotted in Fig. (4.9c) for horizontal cylinder of 16 mm 
diameter. Even though the results obtained in the present work 
are higher than the experimental data they are in better agree- 
ment with experimental data when compared to those obtained 
by Nishikawa et.al. (1976). More importantly, in the present' 
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analysis, the wall is assumed as a perfect black surface 
which can hardly be . met with real material surfaces. 

Nishikawa. et.al. (1976) did not mention the emissivity of 
the cylinder surface. The assumption of emissivit 3 ^ equal to 
1.0 might be one of the reasons for the over prediction 
of average heat transfer coefficient when compared to the 
experimental data. Generally the emissivities of oxidized 
metallic surfaces are in the range of 0.8 to 0.9. 

4.6 Suggestions for Future Work s: 

The results obtained in the present work are for water 
at 1 bar pressure. The program listed in Appendix G can be used 
to obtain results for any liquid at any pressure without any 
difficulty, if the property variations can be obtained as 
functions of temperature and pressure in the polynomial form. 

Also, the analysis done in the present work is for 
laminar film boiling. The more commonly encountered turbulent 
film boiling can be studied with the same program, if good 
data is obtained on the spatial variation of eddy diffusivity 
in the vapour and liquid regions. 

The analysis made in the present work is for natural 
convection film boiling and it is yet to be seen and worth 
taring , if the same program can yield results for a mixed 
convection film boiling problem. This can be achieved numeri- 
cally, by imposing a proper velocity (so that it should neither 
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be biased "towards forced convection nor natural convection) 
at the far end in the liquid region. 

In the present analysis, the vapour lying between the 
wall and interface is assumed to be non-participating gas 
for radiation. Eventhough Sparrow (1964) showed that the 
effect of the participating nature of water vapour is totally 
negligible at low pressures, it is worth doing the analysis 
again with more reasonable assumptions (because, Sparrow 
assumed the vapour to be a gray gas). 
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APPENDIX A 


If the original governing equations are in x-y (Cartesian) 
coordinate system and if, for convenience, the x— y coordinate 
system is to he transformed to £(x, y) — n{x, y) coordinate 
system, then accordingly all the partial derivatives should 
also be transformed into the new £ - n coordinate system using 
the following formulae. 

If the original coordinate system is x - y and the new 
coordinate system is £ - n , where g and n are functions 
of x and y and w is any dependent variable, then 


9 W 

3 W 

LL . 

9 W 

9 n 

3 X 

~ 3 £ ’ 

9 X * + 

9 n 

* 3X 

3 W 

_ 9w 

i* + 

9 W 

9 n 

3 y 

9 £ * 

9y 

9 ri 

* 9y 

jsfw 

9X 2 

j. 2 w 
9 £ 2 

(LI) 2 
v 3x' 

+ 2 

,2 

3 W 

9 £ 4 9 n 


2 2 
3 % 3 a) 3 n 

• o *** ' o 

3 £ 3X Z 3 H 3X^ 


2 

3 £ 9 n 9 w 

3X • ax + an 2 * 


(Al) 

(A2) 


(A3) 


a 2 

3 W 


9 £ 

i i + 

2 

3 w 

r 9 5 

*■ 3 X * 

9 n 

_ 9 _£ 9 _n 

3 X 3 y 

= 35 2 ' 

3 X * 

9 y 

9 £ 3 n 

9y + 

sy * ax 


9 2 W 

9JL 

9 n 

. 9W 

. 3 2 

9 W 

2 

9 n 


9n 2 

* 3 X 

* 9 ? 

9 £ * 

3 X ay 

9 n 

* 9 x 3 y 


(A4) 



A-2 


3 2 w 

sy 2 


.? _W 0 3^w 

' rs ir / 


3 C 


2 ay- 


+ 2 

a c 3;n 


a 


*- tffp. ; W ■&, * 

— . i*? ♦'» 

’r.#, 


M £ji + .Lw. 

a 2 
3 n 


ay 


ay 


(§f)‘ 


•■••JLAZ 
2 


7 5 

**'U :*.■<.■ K ■r' w.yJV 
t./ 


(A5) 


Applying these formulae to the present problem. 


Initially it was x— y coordinate system. New coordinate 
system is x - n 
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Other partial derivatives are irrelevent to this 


problem 



APPENDIX B 


Mass Conservation at the Interface ; 

Mass conservation equation in transient state can 
be derived as follows. 

By taking a control volume at the interface as shown 
in Fig. (2.3) we can write the mass balance as follows s 

Mass of vapour and liquid leaving the 
control volume per unit time 

= U v .d6 - V^.dx + V L dx - U L .d6 (Bl) 

Change in the mass of vapour and liquid in 
the control volume per unit time 

~ p l * ft * dx " p v * If * dx 

(l) and (2) should be equal for mass conservation 




Making use of this, "the equation (2.23) which governs the 
change in boundary-layer thickness with time can be derived 
as follows. By taking a control volume as shown in Fig. (2.4) 
we can write the mass balance as follows. 

Rate of vapour generation — Net flow rate of vapour 

= Rate of change in mass in control volume (B5) 
Now, we shall derive the expression for net flow rate of 
vapour in a differential vapour boundary layer element. By 
observing the same fig. (2.4) the net flow rate in x-direction 
over a differential length dx can be written as 


^x+dx ” *x 


where < 


A 


x+dx 


S+d 6 

= ^ p v^v I x+dx 

6-t-d 3 


. dy and m = / p U I dy 
J x j q v v'x J 


m 


x+dx 


- *x - { p v u v l x«bc dy - { P v U v4 dy 

6+d6 a 

= / [p U | +v~(pU).dxj +. . . ]dy 

J L v v'x 3x v v 1 J 


0 


-/put dy 
J v v'x J 


0 

6+d6 


6+d6 


l -k (f v V- dx - d y + / c v u v d y 

U 6 


6+d 6 . 6+d6 

= - / — (p V ) dx.dy + / p U .dy 
J q jy v v v' J J v v J 


( * . * From continuity equation 


3 


a 






i 


APPENDIX C 

All the thermophysical properties are functions of 
temperature. The equations for the variation of properties 
are obtained from Sabhapathy (1980) , Kothandaraman (1975), 
and steam tables - ERA (1967). All the property values are 
taken in SI units. 

The equations that are used for the calculation of 
thermophysical properties are? 

P v = 219. 5/T ' T in °K 

ky. = (80.4 + 0.407 (T-273.0) ) x 10 -7 for T < 973°K 

= (366.7 + 0.523 (T-973)) x 10" 7 for T > 973°K 

K v = 0.0022 + 0.000098 (T-350) for T < 650°K 

= 0.058 + 0.000145 (T-650) for T > 650°K 

C = (0.4432 - 2.8532 x 10" 5 T + 1.9848 x 10~ 7 T 2 

•^v 

- 6.6372 x 10” 11 T 3 ) x 4186.8 

e L » 972.0 + (353.0 - T) x 0.69984 . 

v L = (4.701 + 0.0916 (333 - T)) x 10~ 4 for T < 333°K 
= (3.556 + 0.05725(353 - T)) x 10“ 4 

for 333°K < T < 353°K 
= (2.821 + 0.03518 (373.0 - T)) 10~ 4 for 

for 353°K < T < 373 °K 
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K l = (668.7 - 0.87 (353-T) )10“ 4 for 333°K < T < 353°K 

= (628.0 - 1.16 (333 - T)) 10~ 4 for T < 333°K 

= (680.4 - 0.555 (373. 0-T)) 10~ 4 for 353°K<T<373°K 

C p = (2.1397 - 9.6814 x 10~ 3 x T 

+ 2.6854 x 10~ 5 x T 2 - 2.4134 x 10” 8 x T 3 ) x 4186.8 

The density is in kg/m 3 , coefficient of viscosity in kg/m-s, 
the thermal conductivity in w/m°K and the specific heat in 
j/kg°K. For non-dimensional form of these property variations 
refer to the program listing. 
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Here, an approximate stabilit3r criterion which limits 
the allowable timestep is given. Eventhough, the criterion 
given below applies only to single-phase flows, it is inte- 
resting to note that it predicts the timestep satisfactorily 
to the present two-phase flow problem. 

The stability requirement of a two-dimensional ad vect ion- 
diffusion equation with upwind differencing was given (Patrick 
J. Roache, 1972) as 


At < 


2a ( U_ „ _1_ ) + M + M 

v 2 2 ' Ax Ay 

Ax Ay y 


Substituting various values in the above inequality (ax = 0.001 m s 

Ay = 5.45 x 10“ 5 m, a = 0.0692/3600, |u| = 0.5661 x 20 m/sec, 

V s 

1 V | = 0.005832 x 20 m/sec.) we can write it as 


At < 


~ 38.44 + 12944.0 + 11320.0 + 2139.0 

< 2645174 

< 37.8 micro seconds 


And the allowable timestep found numerically for the set 
of values given above is 20 micro seconds. So, the predicted 
value from the above inequality is comparable with the 
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actual timestep stipulated by stability of the scheme. By 
looking at the order of magnitude of each term in the above 
inequality it can be concluded .-that Ay has a pronounced 
effect (two-fold) on At.' So, ap we. proceed further seeking 
solution at downstream regions the dimensional value of Ay 
is increasing and this could lead to increase the, allowable 
timestep at downstream regions. This is consistent with the 
remark made in section (3.3). Eventhough at downstream regions, 
U is larger (thus decreases the allowable timestep at down- 
stream regions) its effect is offset bjr the decrease in V 
and a dominant two— fold effect of increase in Ay at downstream 
regions. The net effect is to increase At at downstream 
regions. 

It was mentioned that the timestep was limited by 
stability on vapour-side. This is easy to understand because 
the velocities in vapour-region are much larger than the 
velocities in liquid region. Thus the last two terms in the 
inequality are small in liquid region which implies that a 
relatively large time step can be used in liquid region. 



APPENDIX E 


Transformation of non-linear abscissa (used by 
Nishikawa et.al. (1976)) to a more convenient y/6 coordinate 
is shown in this appendix. Nishikawa 1 et.al. used the non- 

/v 

linear abscissa which is defined as follows. 
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and T, the temperature at any location y, is assumed to be 
linear function of y, for brevity (See Fig. (4.2b) for the 
actual temperature variation) . 
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APPENDIX F 


Number;, 


and 


Nu 

average 

h = h 


h L — i r L 

~ where h = / h dx 

vs 0 


convection + h radiation ~ h c + h r 


1/4 

Also h c x is constant irrespective of x, because of the 

similarity shown by Nishikawa et.al. (1976) for no radiation 

film boiling problem, h^x^^ is calculated for various x 

values in the present work and it is found to be fairly 

constant away from the leading edge. The deviation of hx 1//Zf 

at the leading edge from the constant value at the downstream 

is not more than 6 percent. This is because of the leading 

edge singularity problems in most of the numerical schemes. 

However, this can be totally eliminated by reducing the 

grid size near the leading edge and thus dividing the length 

of the plate into more number of divisions. In the present 
1/4 

work h x ' is approximated to be constant all along the 
length of the plate and thus the local conductive-convective 

fici 
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heat transfer coefficient h c is calculated as follows 
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hc x 

” 7 ^ 

Constant (a value taken at the end of the plate] 

x 1 ' 

(FI) 
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Also, it was observed that the value of h x 1 / 4 is constant 

o 

within the accuracy of 8 percent even when radiation is 
included. Thus, even when radiation is included in the 
numerical scheme, h^ is calculated from the equation (Fl)» 
Hereafter, we refer to h^ as the local conductive—convectivc 
heat transfer coefficient when radiation is included and h 

co 

as the local conductive-convective heat transfer coefficient 
when radiation is absent. 


2. Comparison of Results with those of Sparrow 

Sparrow had given the ratio of q/q Q as a function of X = 

[ g pp X /(16vK T w - . x for various values of 

]' 


and 


N-, 


h 16 v K (T - T ) 1/3 

. r r w vs - ' 

K L * 

g P L X 

o 2 

h = cr(T^ + T)(T+T) 
r w vs' v w vs' 

0.84 C (T - I ) 

X [ 1 + — ^ ] 

X P„ 


For T - T = 800 K, N, = 0.31696 
w vs 1 

(after substituting values of all other parameters at satura- 
tion temperature in the expression for N^) 

and X * 200 corresponds to x = 1,1 cms. 

For X = 200 and 21 °*31 q/q Q was given by Sparrow as ~ 1-92 
From the present analysis . q/q Q = ( h c + ^r^/^co — 1*66, 



3 * &22&2L LB e ^orAccou^bin^ Radiation : 

Bromley’s thumb rule for accounting radiation in 
saturated boiling can be derived as follows. 


* 


m a h. AT 

therefore, m*/A ~ h/h 


and A a’h aT 
co 


CO 


Also, xh a 5^, 

therefore, A*/* ~ <5 /£> ~ h/h 

1/3 


{ / { *= (h co /h) 


'co 


and 


h 


But 


K/6 


h + h = JL 
C r 6* 


CO 




-r * h 

r 


, + >V - h„„ + 

0 

This implies that, 

h = h + h = h ( -——£2—. ^1/3 

r> r» ^ \ 


therefore, h = h . — v + h = h (— + h 

co s * V co v h ; + n r 


h. 


h. 


co v h + h 
c r 


+ h 


h c — ■£££— 'j 1/3 

co ^ h c + h r ; 



= , h co \l/3 ■ 

h ^ h + h ' 

co c r 


(F2) 


This is the condition to be verified to affirm Bromley s 
intuitive rule for accounting radiation. Both the LHS and 
RHS quantities in the equation (F2) are 'evaluated, using the 
results of the present work, and compared. Comparison of 
these two quantities is shown in table. 



Table F-l; Verification of Bromley s Rule 
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4. Transformation of the Prpqpn+ 

Horizontal Cylinder^ — " — — 


h oy Under ~ h c cylinder + had 


1 c cylinder = h / • 

c flat plate (sin 

= h c (Sin O) 1 / 4 

r _ 2 n r D/ S 

c cylinder ~ tcD f h c cylinder 

? ^°/2 


(sin e) 444 


dx V, 


? 71 D / 2 7 

TtD 4 ^ c (Sin ©)■*■ „dx 


1 

'c cylinder nD 


tcD/2 ,, 

/ h c (sin (x/R)) 1 / 4 



g sin © 


X hx lAf /2 r 

TED * V J I- 


il /4 


^ dx 


h cylinder • 1)1/4 = ^ . D 1 /'* (h x 1 / 4 ) { 


Now define, 


/ ^ ] V 4 #dx j 

2x/nD = x s x 1//4 = x 1 / 4 (-) 1 / 4 


h D l/4 - (h v- 1 / 4 ^ n l/4 , 2 xl/4 

c cylinder ° u ~ x ' D * (^Er 


Sin tcx 


f 

0 x 


]l/4 


- h c x 1/4 (f) 1/4 } Sinh^ .ta 

0 x' 

- _L_ h x iA 

" 1.03 ,n c 
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Cylinder * D cylinder * ] ° lA + K» d *^ 4 


nsj *0 * 1/4 + h , d 1 / 4 


(F3) 


The following table shows the error in the prediction of 

average heat transfer coefficient for cylinders of various 
diameters. 




Table F-2s Comparison of present results with experimental data 
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Fig. 2. 3 Control volume at the interface 
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Fig. 2.4 Control volume in vapour region 
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Fig. 4.4 (c ,d ) Mass conservation at the interface 
Fig. 4.4 (e) Mass conservation for a flat interface 

Fig. 4 4 (f) Sign of V velocity near the interface 
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Fig. 4.6 Individual effect of properties of vapour on 
NU Aver aqe * n saturated film boiling 
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